The effect of disorder on the critical temperature of a dilute hard sphere gas 
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We have performed Path Integral Monte Carlo (PIMC) calculations to determine the effect of 
quenched disorder on the superfluid density of a dilute 3D hard sphere gas. The disorder was 
introduced by locating set of hard cylinders randomly inside the simulation cell. Our results indicate 
that the disorder leaves the superfluid critical temperature basically unchanged. Comparison to 
experiments of helium in Vycor is made. 
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Recently there have been several theoretical studies of 
the superfluid transition (Bose condensation) of a sys- 
tem of homogeneous hard spheres. For certain properties 
the model of hard spheres can describe the Bose-Einstein 
condensates of alkali metal atoms and liquid 4 He. The 
ratio of the superfluid transition temperature to that of 
an ideal gas T c /T c q was determined with numerical cal- 
culations jjj to be larger than unity over a wide range 
of densities. The enhancement has been confirmed by 
several other studies H|| however the magnitude of the 
enhancement as a function of density is still controversial. 
It has also been proposed that T c can be enhanced with 
respect to an ideal gas by the introduction of quenched 
(i.e. static) disorder Experimental results for 4 He in 
Vycor show that, for some concentrations, the enhance- 
ment is as large as 300 %. Here, we study numerically 
how a low density hard spheres gas is affected by ran- 
dom geometry. Although there have been many calcu- 
lations for disordered lattice models ||,||, to our knowl- 
edge, there are no other calculations tackling this prob- 
lem in the continuum. 

Both helium gas and hard spheres system fall into the 
universality class in which the critical exponents of the 
superfluid fraction, ip, and the correlation length, v, sat- 
isfy cp = -v = -0.67. In this class, the critical exponent 
of the specific heat, a is negative and according to the 
Harris criterium j7j], the presence of weak uncorrelated 
disorder will not change the exponents. To examine 
both this change in v and the possible variation in T c 
with the presence of quenched disorder, one can observe 
the behaviour of helium when absorbed in a network of 
porous silica |8| |l^] . This material is produced with var- 
ious porosities. For example, in aerogel the volume of 
the voids in the system is around 90-95%, in xerogel or 
porous gold it is 60%, and in Vycor the porosity is ~ 
30%. Experimentally, it is found that the introduction of 
disorder decreases both the critical temperature and the 
fraction of the helium atoms that decouple from the os- 
cillator measurement of the superfluid. In addition, the 
critical exponents can be changed: for samples with large 
porosity, v has a value greater than of bulk helium case 
0,0. However, in Vycor |1(J and porous gold jL4|, the 



exponents are the same as in bulk 4 He. 

Since a system of homogeneous hard spheres belongs to 
the same universality class as bulk liquid 4 He, we tried to 
understand the changes in the critical temperature and 
critical exponents in silica gels by performing Path Inte- 
gral Monte Carlo (PIMC) calculations jl5| for a system 
of hard spheres with quenched disorder. A simple model 
of disorder was constructed by placing, at random, N c 
hard cylinders parallel to the edges of the cubic periodic 
simulation cell. We used hard cylinders rather than a 
more realistic model of Vycor or aerogel, because such 
a model is much easier to treat computationally, in part 
because the absence of an attractive interaction means 
that there will not be a "dead layer" of 4 He next to the 
substrate which would impose a computational burden 
without participating in superflow. We chose the posi- 
tions of the cylinders by doing a simulation of a 2D hard 
quantum hard disks. One snapshot was used to place the 
cylinders along the x-axis, another independent snapshot 
for the y-axis and a third for the z-axis. Thus there were 
always equal numbers of cylinders along each axis. 

The hard spheres are chosen to have diameter a appro- 
priate to 4 He (2.2 A) so as to make comparison straight- 
forward, though we will also use dimensionless units. We 
demand that the hard spheres remain a distance d away 
from the center of each of the cylinders. The diameters 
of the cylinders were varied to study the effect of their 
number and/or of the total excluded volume. We con- 
sidered two densities: pa 3 = 0.0107 and 1.07 x 10~ 5 . 
Here, p = N s /L 3 , N s is the number of hard spheres and 
L the length of the simulation cell. A summary of the 
conditions used is given in Table I. Each simulation was 
repeated for five different realizations of the disorder. 

There are four dimensionless parameters characteriz- 
ing this system in the thermodynamic limit: the density 
of hard spheres pa 3 , the fraction of excluded volume, an 
upper bound is p = ■nN c (d/L) 2 (an upper bound be- 
cause if two cylinders are close to each other, some of the 
excluded volume is counted twice); the ratio d/a; and 
the ratio of the temperature to that of free boson tran- 
sition temperature at the same density (T c o=0.401 K at 
pa 3 = 0.0107, not taking into account excluded volume.) 



First, let us examine how the precise geometry of the 
disorder affects the superfluid fraction. Fig. 1 displays 
p s /p, for system II, as a function of T/T c0 . The super- 
fluid density pl| is computed by calculating the mean 
squared winding of the paths around the periodic bound- 
ary conditions. The crosses indicate the results for each 
disorder at a given temperature, and the lines are the 
smoothed averages of the superfluid fractions for the five 
cases. All the values of p s /p are within a standard de- 
viation of the average values, indicating that the super- 
fluid fraction does not depend much on the particular 
geometry chosen and that one can speak of a single T c . 
Even though our systems small they seem to be "self- 
averaging." The behaviour of system II is typical of the 
other three geometries. 
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FIG. 1. The upper curve shows the superfluid fraction for 
system II with N c =9 and iV s =118 The symbols correspond 
to individual realizations, and the curve with error bars to 
the disorder average. The error bars indicate one standard 
deviation. The lower curve shows results for an ideal gas. 

The data in the lower part of the figure were obtained 
by "switching off" the interactions between the hard 
spheres and keeping constant the rest of the simulation 
conditions, to gauge the importance of the interaction on 
the results. We refer to this as the "disordered" ideal gas. 
We did this for systems I, II and III. It can be observed 
that when the hard core interaction is turned off, the 
fraction of molecules belonging to the superfluid at low 
temperatures decreases by a factor of two. In the ther- 
modynamic limit, the effect is even greater; the disorder 
suppresses the superfluidity of the ideal gas. This can be 
easily understood since, without interaction, there exists 
a lowest energy state, localized in the cavity with the 
largest volume, that will be occupied by all the bosons. 
Thus, when N s —> oo, one will have BEC but not super- 
fluidity. In our simulations, we observed many exchanges, 
but in the absence of interaction, the exchange loops re- 



mained in a single pocket. The winding numbers were 
much less than in the interacting case and were consistent 
with scaling to a non-superfluid state at all temperatures. 

Another effect we observe in our simulation is the re- 
duction of superfluid density with increasing disorder. 
That can be seen in Fig. 2 for the different systems 
with N s ~ 60. The tortuosity, which we define as the 
fraction of atoms participating in superflow at zero tem- 
perature, x = Ps(T = 0)/p given in Table I, is found to 
be independent of N s and 0.51 < \ < 0.74. These values 
are reasonably consistent with those measured in Vycor 
X = 0.76±0.01 ||. The greater the excluded volume (sys- 
tem III vs. system II) for the same type of cylinders, the 
lower x- For the same excluded volume (systems I, II and 
IV), the superfluid fraction is greater when the number 
of cylinders is smaller, if the diameter of the hard spheres 
is comparable to that of the cylinders (d ~ a) . 
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FIG. 2. The superfluid density for the pure system (pluses) , 
system II (crosses), I (stars), IV (open squares) and III (full 
squares). Here 54 < N„ < 64. The same ordering of the 
superfluid density occurred for other values of N s . The error 
bars are shown for one of the systems, and are comparable in 
all other cases, except in the clean system, which are a factor 
of two smaller. 

To obtain an estimate of the superfluid transition tem- 
perature for each system, we perform a finite size scaling 
analysis by fitting N^ 3 p s /p with fl6|| : 

A + BN ( S 1/3 ^ (T - T c ) + CN~ A/{3 "\ (1) 

This expression assumes the exponents for the superfluid 
density and the correlation length are related by (tp = 
-v). The fit determines the critical temperature, T c , and 
the exponent for the correlation length, v. We assumed 
that A = 1/2, but its precise value has very little effect 
on the results. The range of reduced temperatures Ti < 
T < T 2 ( Ti and T 2 given in Table II), included in the 
fit was adjusted to ensure that the superfluid density was 



linear. The fit for system I is given in Fig. 3. Shown in 
Table II are the transition temperatures obtained with 
two different fits: one with the i/s corresponding to the 
universality class of the hard spheres (helium) , and other 
choosing v to minimize x 2 . Note that T c is independent 
of the value of v chosen. The introduction of the cylinders 
does not change T c much from an homogeneous ideal 
gas at the same density. Recall that without disorder 
T c /T c0 = 1.057 ± 0.003 § (see figure 4). 
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FIG. 3. Scaled superfluid density for system I. The symbols 
refer to N s = 19, 54, 99 and 153. The dotted lines are the fits 
and their intersection marks the critical temperature for this 
system, in this case T c /T c o= 1.02 ± 0.05. The error bars are 
given for iV,,=153, and are at larger at least by a factor of two 
than in the other cases. 

As the density of cylinders increases (system II vs. I) 
or their diameter increases (system I vs. IV), the criti- 
cal temperature increases. This is most likely because as 
the available space for the hard spheres decreases, their 
effective density increases. On the other hand, for the 
same kind of cylinders (systems II and III), T c does not 
change when the excluded volume does. Figure (4) shows 
the comparison with the bulk system. We see that within 
error bars, disorder does not increase the transition tem- 
perature, but on the other hand, it is not reduced either. 

With respect to the universality class of the systems 
under consideration, our conclusions are severely limited 
by the statistical quality of the parameters obtained from 
the fit to Eq. (1). For the most disordered system (III), 
the one with a excluded volume of 0.33, we cannot draw 
any conclusion about the value of v, so we cannot say 
if the universality class has changed. Luckily, this fact 
is irrelevant for the critical temperatures. In the other 
three systems, the situation is better. The error bars of 
the critical exponent for the correlation length are small 
enough to conclude that systems I and II belong to the 
same universality class as pure helium (for which v = 
0.67). On the other hand, in system IV v = 1.32 ± 0.08. 



This value is clearly different of 0.67, a result reinforced 
by the dramatic decrease in % 2 when we go from v = 0.67 
to the best fit of v — 1.32 (sec table II). This system might 
be in a different universality class from pure helium. Such 
an increase of v has been experimentally observed in he- 
lium absorbed in aerogel and xerogel ]l(|[il"|]i"3| ] . 

We now compare these calculations to the measure- 
ments of the transition temperature in Vycor. There are 
several important differences. In Vycor the first helium 
atoms added to the system are localized on the surface of 
the silica and do not contribute to a superfluid response 
(the dead layer). The next helium atoms form a surface 
layer around the silica strands which is locally 2D but 
connects up at a longer length scale as a 3D network. 
The pores are estimated to be in the range of 40 to 80 
A. For the sample of ref. ||J^] we estimate the excluded 
volume, including that of the dead layer, is [i = 0.772. 
This is a much larger excluded volume factor and a much 
larger pore size than in our simulations. The tortuosity 
of Vycor is estimated to have a value of x = 0.76 ±0.01 
similar to that of our systems. 
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FIG. 4. The superfluid transition temperature versus the 
density. The solid line is the calculations of ref. [1]. The 
o are the data shown in fig. 3 of ref [4]. They have been 
corrected for the free volume, tortuosity, and effective mass. 
The solid symbols with error bars are the present calculations. 
The dashed line is for bulk liquid 4 He. 

Figure 4 shows the comparison of the current results 
with the experiment. To make this comparison, we cor- 
rected the densities for the excluded volume in both sets 
of data. Plotted on the horizontal axis is the density of 
helium in units of hard sphere radii, defined as the num- 
ber of free atoms (not including the dead layer) divided 
by the free pore volume. The vertical axis is the tran- 
sition temperature divided by the ideal transition tern- 



perature also taking into account the corrected density. 
The upper curve is the estimate for pure hard spheres of 
rcf 0. The transition temperature in Vycor is reduced 
from the bulk hard spheres value much more than in our 
calculation, by up to a factor of three, most likely be- 
cause of the density variation within a pore and between 
pores. Assuming a pore size of 40A, for densities such 
that pa 3 /(l — /j,) > 1.1 x 10~ 3 , one has more than one 
free atom/per pore, thus for these densities, bose con- 
densation can take place inside a single pore. However, 
not until the pores are nearly full, is our model of strictly 
repulsive disorder relevant for Vycor because of layering 
within a pore. This effect, as well as the significant dif- 
ference in the ratio of pore size to hard sphere diameter, 
likely accounts for the different reductions in the tran- 
sition temperature at low density. Also plotted on fig. 
(4) is the transition temperature of bulk 4 He versus den- 
sity, showing essentially the same reduction as for 4 He in 
Vycor at the same density. 

In conclusion we have performed PIMC simulation to 
determine the superfluid density for hard sphere system 
outside of randomly placed cylinders. We have an over- 
all reduction in the superfluid density, however the tak- 
ing into account the reduced free volume, the transition 
temperature is reduced only slightly if at all. 
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System 




A* 


d/a 


P,(T = 0)/p 


N c 


N s 


I 


1.07 x 10 -2 


0.128 


1 


0.68±0.10 


6 


19 












12 


54 












18 


99 












24 


153 


II 


1.07 x 10 -2 


0.128 


1.5 


0.74 ± 0.05 


3 


23 












6 


64 












9 


118 












12 


182 


III 


1.07 x 10" 2 


0.33 


1.5 


0.51±0.19 


9 


29 












15 


62 












21 


102 












30 


174 


IV 


1.07 x 10" b 


0.128 


15 


0.56±0.09 
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23 
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64 
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118 












12 


182 



System 


Ti/T c o 


T 2 /T c0 


T c /T c o 


u 


x 2 


I 


0.72 


1.12 


1.02 ± 0.05 


0.67 (fixed) 


0.8 


I 


0.72 


1.12 


1.02 ± 0.05 


0.66 ± 0.05 


0.8 


II 


0.72 


1.20 


1.15 ± 0.05 


0.67 (fixed) 


1.4 


II 


0.72 


1.20 


1.15 ± 0.05 


0.71 ± 0.05 


1.4 


III 


0.72 


1.25 


1.15 ± 0.20 


0.67 (fixed) 


0.7 


III 


0.72 


1.25 


1.12 ± 0.75 


1.30 ± 0.5 


0.4 


IV 


0.62 


1.25 


1.17 ± 0.05 


0.67 (fixed) 


14.5 


IV 


0.62 


1.25 


1.17 ± 0.17 


1.32 ± 0.08 


1.6 



TABLE II. Transition temperature T c and exponent v as 
determined by the fit to Eq. (1). Ti < T <T2 is the fit range, 
\ 2 is the reduced fit residual. The number of data points in 
each fit was 32. 
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TABLE I. Different cylinder arrangements and simulation 
conditions used in this work. Systems II and IV have the 
same randomness (placement of the cylinders) 



